function smp = srf_DensityMap(hfile, opts)
% SRF::DensityMap  - generate a vertex density map
%
% FORMAT:       dsmp = srf.DensityMap([opts])
%
% Input fields:
%
%       opts        optional settings
%        .nei       use distance of n-th. neighbor (default 1)
%
% Output fields:
%
%       dsmp        density SMP (in average mm distance to neighbors)

% Version:  v0.7b
% Build:    7092614
% Date:     Sep-26 2007, 2:55 PM CEST
% Author:   Jochen Weber, Brain Innovation, B.V., Maastricht, NL
% URL/Info: http://wiki.brainvoyager.com/BVQXtools

% argument check
if nargin < 1 || ...
    numel(hfile) ~= 1 || ...
   ~isBVQXfile(hfile, 'srf')
    error( ...
        'BVQXfile:BadArgument', ...
        'Invalid call to ''%s''.', ...
        mfilename ...
    );
end
if nargin < 2 || ...
    numel(opts) ~= 1 || ...
   ~isstruct(opts)
    opts = struct;
end
if ~isfield(opts, 'nei') || ...
    numel(opts.nei) ~= 1 || ...
   ~isa(opts.nei, 'double') || ...
   ~any((1:5) == opts.nei)
    opts.nei = 1;
end

% get content
sc = bvqxfile_getscont(hfile.L);
bc = sc.C;
c = bc.VertexCoordinate;
nc = size(c, 1);
n = bc.Neighbors;
n2 = n(:, 2);

% for other than first neighbor
if opts.nei > 1
end

% get average distance to neighbors
ad = single(zeros(nc, 1));
md = single(zeros(nc, 1));
sd = single(zeros(nc, 1));
for vc = 1:nc
    vd = sqrt(sum( ...
        (repmat(c(vc, :), [numel(n2{vc}), 1]) - c(n2{vc}, :)) .^ 2, 2));
    ad(vc) = mean(vd);
    md(vc) = min(vd);
    if numel(n2{vc}) > 1
        sd(vc) = std(vd);
    end
end

% create smps
smp = bless(BVQXfile('new:smp'), 1);
smpc = bvqxfile_getcont(smp.L);
smpc.NrOfVertices = nc;
smpc.NameOfOriginalSRF = sc.F;
smpc.Map.LowerThreshold = double(min(ad));
smpc.Map.UpperThreshold = double(max(ad));
smpc.Map.DF1 = round(3 * (size(bc.TriangleVertex, 1) / nc));
smpc.Map.Name = 'Average distance to %d. neighbors';
smpc.Map.SMPData = ad;
ad = 1 ./ (ad + eps);
smpc.Map(2) = smpc.Map(1);
smpc.Map(2).LowerThreshold = double(min(ad));
smpc.Map(2).UpperThreshold = double(mean(ad) + 2 * std(ad));
smpc.Map(2).Name = 'Vertex density w.r.t. %d. neighbors';
smpc.Map(2).SMPData = ad;
smpc.Map(3) = smpc.Map(1);
smpc.Map(3).LowerThreshold = double(min(md));
smpc.Map(3).UpperThreshold = double(mean(md) + 2 * std(md));
smpc.Map(3).Name = 'Minimum distance to %d. neighbors';
smpc.Map(3).SMPData = md;
smpc.Map(4) = smpc.Map(1);
smpc.Map(4).LowerThreshold = double(min(sd));
smpc.Map(4).UpperThreshold = double(mean(sd) + 2 * std(sd));
smpc.Map(4).Name = 'Local distortion (std) in distance to neighbors';
smpc.Map(4).SMPData = sd;
smpc.NrOfMaps = 4;
bvqxfile_setcont(smp.L, smpc);
